##Create .log file
replog <- file("SPL_ReplicationCode_Log.log") 
sink(replog, append = TRUE, type = "output")
sink(replog, append = TRUE, type = "message")
cat(readChar(rstudioapi::getSourceEditorContext()$path, 
             file.info(rstudioapi::getSourceEditorContext()$path)$size))

########################################
########################################
########################################
########################################
##Working directory should be set to location of this R file.
########################################
########################################
########################################
########################################

########################################
########################################
########################################
########################################
##Replication code for statistical analysis reported in "Complementary Partners? Attitudes toward Multi-Actor Development Projects in the Democratic Republic of Congo"
##Austin Strange, Elizabeth Plantan, and Wendy Leutert
##Last updated: 8 August 2023
########################################
########################################
########################################
########################################
 
##Load packages
library(arsenal)
library(broom)
library(cobalt)
library(dplyr)
library(ggplot2)
library(MASS)
library(plyr)
library(stargazer)
library(tidyr)
library(vtable)
library(xtable)

########################################
########################################
########################################
########################################
##Load, merge, and clean data
########################################
########################################
########################################
######################################## 

##Read cleaned data containing variables used in manuscript analysis (including manuscript and supporting materials)

##Full dataset used for primary analysis
d <- read.csv("SPL_4August23.csv")

##Update treatment groupings into factors with appropriate levels
d$TreatGroup <- as.factor(d$TreatGroup)
d$TreatGroup <- relevel(d$TreatGroup, ref = "CooperationControl")
d$brand_pooled <- as.factor(d$brand_pooled)
d$brand_pooled <- relevel(d$brand_pooled, ref = "Coop")
d$ext_pooled <- as.factor(d$ext_pooled)
d$ext_pooled <- relevel(d$ext_pooled, ref = "NoExt")

##Subsetted dataset of respondents in cooperation (multi-actor) treatment groups
dext <- read.csv("SPL_Cooperation_4August23.csv")

##Subsetted dataset of respondents in cooperation (multi-actor) treatment group with negative externality
dextneg <- read.csv("SPL_CooperationNeg_4August23.csv")

##Subsetted dataset of respondents in treatments group with no externalities
d.noext <- read.csv("SPL_NoExt_4August23.csv")

##Update treatment groupings into factors with appropriate levels
d.noext$TreatGroup <- as.factor(d.noext$TreatGroup)
d.noext$TreatGroup <- relevel(d.noext$TreatGroup, ref = "CooperationControl")

########################################
########################################
########################################
########################################
##Manuscript figures and tables 
########################################
########################################
########################################
########################################

########################################
########################################
##Table 1: Satisfaction levels and difference-in-means 
########################################
########################################

#Compare differences in means by treatment group
t1a <- tidy(t.test(d$PosTre6_Brand_Satisfn[d$TreatGroup == "CooperationControl"], d$PosTre6_Brand_Satisfn[d$TreatGroup == "INGOControl"], data = d))
t1b <- tidy(t.test(d$PosTre6_Brand_Satisfn[d$TreatGroup == "CooperationControl"], d$PosTre6_Brand_Satisfn[d$TreatGroup == "ChinaControl"], data = d))
t1c <- tidy(t.test(d$PosTre6_Brand_Satisfn[d$TreatGroup == "CooperationPositive"], d$PosTre6_Brand_Satisfn[d$TreatGroup == "INGOPositive"], data = d))
t1d <- tidy(t.test(d$PosTre6_Brand_Satisfn[d$TreatGroup == "CooperationPositive"], d$PosTre6_Brand_Satisfn[d$TreatGroup == "ChinaPositive"], data = d))
t1e <- tidy(t.test(d$PosTre6_Brand_Satisfn[d$TreatGroup == "CooperationNegative"], d$PosTre6_Brand_Satisfn[d$TreatGroup == "INGONegative"], data = d))
t1f <- tidy(t.test(d$PosTre6_Brand_Satisfn[d$TreatGroup == "CooperationNegative"], d$PosTre6_Brand_Satisfn[d$TreatGroup == "ChinaNegative"], data = d))

#Store outputs needed to produce Table 1, and produce equivalent table
#Treatment group sizes
ns <- d %>%                           
  group_by(TreatGroup) %>% 
  dplyr::summarize(count = n())

#Mean satisfaction levels
ms <- aggregate(d$PosTre6_Brand_Satisfn, list(d$TreatGroup), FUN=mean) 

#Differences in means and p-values
t1 <- cbind(ns,ms); t1 <- t1[,c(1:2,4)]; colnames(t1)[2:3] <- c("N","Satisfaction Level"); t1$DiffMeans <- ""; t1$pvalue <- ""

#Remove baseline group not used for this hypothesis testing; rearrange treatment groups in line w/ Table 1
t1 <- t1[t1$TreatGroup != "Control_Control",]; t1 <- t1[c(1,7,2,6,9,4,5,8,3),]

#Incorporate differences in means and p-values
t1$DiffMeans[c(2:3,5:6,8:9)] <- c(t1a$estimate,t1b$estimate,t1c$estimate,t1d$estimate,t1e$estimate,t1f$estimate)
t1$pvalue[c(2:3,5:6,8:9)] <- c(t1a$p.value,t1b$p.value,t1c$p.value,t1d$p.value,t1e$p.value,t1f$p.value)

#Print equivalent of Table 1
print(xtable(t1))

########################################
########################################
##Table 2: Attribution and difference-in-means 
########################################
########################################

#Compare differences in means for attribution of blame (in negative externality scenario) to INGO and Chinese SOE
t2a <- tidy(t.test(dextneg$PosTre6_Ext_RespPrim_INGO, dextneg$PosTre6_Ext_RespPrim_ChinaSOE, data = dextneg))

#Compare differences in means for attribution of responsibility (in negative externality scenario) to INGO and Chinese SOE
t2b <- tidy(t.test(dext$PosTre6_Ext_RespSol_INGO[dext$ext_pooled == "Neg"], dext$PosTre6_Ext_RespSol_ChinaSOE[dext$ext_pooled == "Neg"], data = dext)) 

#Compare differences in means for attribution of responsibility (in positive externality scenario) to INGO and Chinese SOE
t2c <- tidy(t.test(dext$PosTre6_Ext_RespScs_INGO[dext$ext_pooled == "Pos"], dext$PosTre6_Ext_RespScs_ChinaSOE[dext$ext_pooled == "Pos"], data = dext))

#Store outputs needed to produce Table 2, and produce equivalent table
#Mean attribution 
t2 <- rbind(t2a,t2b,t2c); t2 <- t2[,c(2,3,1,5)]; colnames(t2) <- c("INGO","Chinese SOE","DiffMeans","pvalue")

#Print equivalent of Table 2
print(xtable(t2))

########################################
########################################
##Table 3: Perceived positive impact/effectiveness level difference-in-means 
########################################
########################################

#Compare differences in means by treatment group, evaluation of INGO
t3a <- tidy(t.test(dext$PosTre6_Brand_Eff_EuroNGO[dext$TreatGroup == "CooperationPositive"], dext$PosTre6_Brand_Eff_EuroNGO[dext$TreatGroup == "CooperationControl"], data = dext)) 
t3b <- tidy(t.test(dext$PosTre6_Brand_Eff_EuroNGO[dext$TreatGroup == "CooperationNegative"],dext$PosTre6_Brand_Eff_EuroNGO[dext$TreatGroup == "CooperationControl"], data = dext)) 

#Compare differences in means by treatment group, evaluation of Chinese SOE
t3c <- tidy(t.test(dext$PosTre6_Brand_Eff_ChinaSOE[dext$TreatGroup == "CooperationPositive"],dext$PosTre6_Brand_Eff_ChinaSOE[dext$TreatGroup == "CooperationControl"], data = dext)) 
t3d <- tidy(t.test(dext$PosTre6_Brand_Eff_ChinaSOE[dext$TreatGroup == "CooperationNegative"], dext$PosTre6_Brand_Eff_ChinaSOE[dext$TreatGroup == "CooperationControl"], data = dext))

#Store outputs needed to produce Table 3, and produce equivalent table
#Mean impact levels
e_ingo <- aggregate(dext$PosTre6_Brand_Eff_EuroNGO, list(dext$TreatGroup), FUN=mean,na.rm=TRUE); e_ingo <- e_ingo[c(1,3,2),] 
e_chn <- aggregate(dext$PosTre6_Brand_Eff_ChinaSOE, list(dext$TreatGroup), FUN=mean,na.rm=TRUE); e_chn <- e_chn[c(1,3,2),]  
t3 <- rbind(e_ingo,e_chn); colnames(t3)[1:2] <- c("Actor/Treatment","Positive Impact"); t3$DiffMeans <- ""; t3$pvalue <- ""

#Incorporate differences in means and p-values
t3$DiffMeans[c(2:3,5:6)] <- c(t3a$estimate,t3b$estimate,t3c$estimate,t3d$estimate)
t3$pvalue[c(2:3,5:6)] <- c(t3a$p.value,t3b$p.value,t3c$p.value,t3d$p.value)

#Print equivalent of Table 3
print(xtable(t3))

########################################
########################################
########################################
########################################
##Supporting information figures and tables 
########################################
########################################
########################################
########################################

########################################
########################################
########################################
##Appendix 3: Descriptive Statistics
########################################
########################################
########################################

########################################
########################################
##Table A1: Summary statistics for regressions
########################################
########################################

#Get statistics for following variables in Table A1
t.a1 <- subset(d, select = c(PosTre6_Brand_Satisfn,sat, BG2_Age, female, Dem3_Edu_if_col, Dem3_Occpn_if_desk, AddOut8_KnowCh, PreT4_Context_PoliKnowle,
                              Knowl9_IMAGINE, act_qrt,PreT4_Eff_INGO,PreT4_Eff_ChinaSOE,AddOut8_WorCovid))
st(t.a1,out='latex')

########################################
########################################
##Table A2: Summary Statistics (Overall and by City)
########################################
########################################

#Get statistics for following variables in Table A2
covsbal <- subset(d, select = c(TreatGroup, BG2_Age,BG2_Sex_1,Dem3_Edu_if_col,Dem3_Occpn_if_desk,Knowl9_IMAGINE,act_qrt,AddOut8_KnowCh,
                                PreT4_Context_PoliKnowle,Dem3_Ownersh_TV,BG2_City))
t.a2 <- tableby(BG2_City ~ ., data = covsbal,test=FALSE)
summary(t.a2, title = "Summary Statistics",text="latex")

########################################
########################################
##Table A3: Balance Table of Core Variables by Treatment Group
########################################
########################################

#Get statistics for following variables in Table A3
t.a3 <- tableby(TreatGroup ~ ., data = covsbal,test=FALSE)
summary(t.a3, title = "Balance Table",text="latex")

########################################
########################################
##Table A4: Pre-treatment Perceptions of Development Actor Effectiveness
########################################
########################################

#Get statistics for following variables in Table A4
pretreat_eff <- subset(d, select = c(PreT4_Eff_WorldBank, PreT4_Eff_DFID, PreT4_Eff_INGO, 
                                     PreT4_Eff_NatiGov, PreT4_Eff_ProvGov, PreT4_Eff_CityGov,
                                     PreT4_Eff_Cong_Pub, PreT4_Eff_EuroComp,PreT4_Eff_ChinaSOE,
                                     PreT4_Eff_CongComp, BG2_City))
t.a4 <- tableby(BG2_City ~ ., data = pretreat_eff)
summary(t.a4, title = "Summary Statistics",text="latex")

########################################
########################################
##Figure A2: Covariate Balance Plots (Figure A3 in conditionally accepted version of manuscript)
########################################
########################################

#Create balance plots for following variables
colnames(covsbal) <- c("Treatment group","Age","Gender","University","Desk job","IMAGINE awareness","IMAGINE neighborhood","Chinese neighbors",
                       "Political knowledge","TV ownership")
covsbal$Gender <- as.character(covsbal$Gender)
covsbal$Female <-ifelse(covsbal$Gender == "Female",1,0)

bp.uni <- bal.plot(covsbal, treat = covsbal$'Treatment group',
                   var.name = "University")
bp.job <- bal.plot(covsbal, treat = covsbal$'Treatment group',
                   var.name = "Desk job")
bp.age <- bal.plot(covsbal, treat = covsbal$'Treatment group',
                   var.name = "Age")
bp.gender <- bal.plot(covsbal, treat = covsbal$'Treatment group',
                      var.name = "Female")
bp.pn <- bal.plot(covsbal, treat = covsbal$'Treatment group',
                  var.name = "Political knowledge")
bp.tv <- bal.plot(covsbal, treat = covsbal$'Treatment group',
                  var.name = "TV ownership")
bp.imag1 <- bal.plot(covsbal, treat = covsbal$'Treatment group',
                     var.name = "IMAGINE awareness")
bp.imag2 <- bal.plot(covsbal, treat = covsbal$'Treatment group',
                     var.name = "IMAGINE neighborhood")

########################################
########################################
##Figures A3 and A4: Reported Reasons for Satisfaction/Dissatisfaction with Project (Figures A6 and A7 in conditionally accepted version of manuscript)
########################################
########################################

#Relevel treatment group factor variable to match figure in manuscript
d.noext$TreatGroup <- relevel(d.noext$TreatGroup, ref = "ChinaControl")

#plot reasons for satisfaction
msat <- d.noext %>%                                       
  group_by(TreatGroup) %>%                         
  summarise_at(vars(sat_qual, sat_health, sat_env, sat_econ, sat_actoreff, sat_integover, sat_citinput),              
               list(name = mean)) 
msat$trt <- msat$TreatGroup; msat <- msat[,c(9,1:8)]
colnames(msat) <- c("Group","TreatGroup","Quality","Health","Environment","Economy","Actor Effectiveness","Integrity","Citizen Input")
msat <- msat %>% gather(Reason, TreatGroup, 'Quality':'Citizen Input'); colnames(msat)[3] <- "Mean"
msat$Mean <- msat$Mean*100
msat$Reason <- factor(msat$Reason, levels=c("Quality", "Health","Environment","Economy","Actor Effectiveness","Integrity","Citizen Input"))
f.a3 <- ggplot(msat, aes(x = Reason, y = Mean, fill = Group)) +
  geom_col(position = "dodge") +
  xlab("Reason") + ylab("Percentage of Respondents") +
  theme(axis.title.x = element_text(size=20),
        axis.text.x  = element_text(size=14, color="black"),
        axis.text.y  = element_text(size=14, color="black"),
        axis.title.y = element_text(size=20),
        legend.title=element_text(size=20),
        legend.text=element_text(size=14))

#plot reasons for dissatisfaction
munsat <- d.noext %>%                                       
  group_by(TreatGroup) %>%                         
  summarise_at(vars(unsat_qual, unsat_health, unsat_env, unsat_econ, unsat_actoreff, unsat_integover, unsat_citinput),              
               list(name = mean)) 
munsat$trt <- munsat$TreatGroup; munsat <- munsat[,c(9,1:8)]
colnames(munsat) <- c("Group","TreatGroup","Quality","Health","Environment","Economy","Actor Effectiveness","Integrity","Citizen Input")
munsat <- munsat %>% gather(Reason, TreatGroup, 'Quality':'Citizen Input'); colnames(munsat)[3] <- "Mean"
munsat$Mean <- munsat$Mean*100
munsat$Reason <- factor(munsat$Reason, levels=c("Quality", "Health","Environment","Economy","Actor Effectiveness","Integrity","Citizen Input"))
f.a4 <- ggplot(munsat, aes(x = Reason, y = Mean, fill = Group)) +
  geom_col(position = "dodge") +
  xlab("Reason") + ylab("Percentage of Respondents") +
  theme(axis.title.x = element_text(size=20),
        axis.text.x  = element_text(size=14, color="black"),
        axis.text.y  = element_text(size=14, color="black"),
        axis.title.y = element_text(size=20),
        legend.title=element_text(size=20),
        legend.text=element_text(size=14))

########################################
########################################
########################################
##Appendix 4: Additional Tests
########################################
########################################
########################################

########################################
########################################
##Table A5: Determinants of Project Satisfaction
########################################
########################################

#Run regressions reported in Table A5
reg_pool_1 <- lm(PosTre6_Brand_Satisfn ~ brand_pooled + ext_pooled + brand_pooled*ext_pooled,data=d); summary(reg_pool_1)
reg_pool_2 <- lm(PosTre6_Brand_Satisfn ~ brand_pooled + ext_pooled + brand_pooled*ext_pooled + BG2_Age + BG2_Sex_1 + Dem3_Edu_if_col + 
                   Dem3_Occpn_if_desk + BG2_City + Knowl9_IMAGINE + act_qrt + AddOut8_KnowCh + PreT4_Context_PoliKnowle +
                   BG2_City,data=d)
stargazer(reg_pool_1, reg_pool_2, single.row = T, digits = 2, star.cutoffs = c(.05, .01, .001,.0001))

########################################
########################################
##Table A6: Determinants of Project Satisfaction (with Neighborhood Fixed Effects)
########################################
########################################

#Run regressions reported in Table A6
reg_pool_1fe <- lm(PosTre6_Brand_Satisfn ~ brand_pooled + ext_pooled + brand_pooled*ext_pooled + + factor(BG2_Qrt),data=d)
reg_pool_2fe <- lm(PosTre6_Brand_Satisfn ~ brand_pooled + ext_pooled + brand_pooled*ext_pooled + BG2_Age + BG2_Sex_1 + Dem3_Edu_if_col + 
                     Dem3_Occpn_if_desk + BG2_City + Knowl9_IMAGINE + act_qrt + AddOut8_KnowCh + PreT4_Context_PoliKnowle +
                   factor(BG2_Qrt),data=d)
stargazer(reg_pool_1fe, reg_pool_2fe, single.row = T, digits = 2, star.cutoffs = c(.05, .01, .001,.0001))

########################################
########################################
##Table A7: Determinants of Project Satisfaction (with Enumerator Fixed Effects)
########################################
########################################

#Run regressions reported in Table A7
reg_pool_e1 <- lm(PosTre6_Brand_Satisfn ~ brand_pooled + ext_pooled + brand_pooled*ext_pooled + factor(enumid),data=d)
reg_pool_e2 <- lm(PosTre6_Brand_Satisfn ~ brand_pooled + ext_pooled + brand_pooled*ext_pooled + BG2_Age + BG2_Sex_1 + Dem3_Edu_if_col + 
                    Dem3_Occpn_if_desk + BG2_City + Knowl9_IMAGINE + act_qrt + AddOut8_KnowCh + PreT4_Context_PoliKnowle +
                    factor(enumid),data=d)
stargazer(reg_pool_e1, reg_pool_e2, single.row = T, digits = 2, star.cutoffs = c(.05, .01, .001,.0001))

########################################
########################################
##Table A8: Support for Multi-Actor: Non-Pooled Data
########################################
########################################

#Relevel treatment group factor variable to match figure in manuscript
d.noext$TreatGroup <- relevel(d.noext$TreatGroup, ref = "CooperationControl")

#Run regressions reported in Table A8
h1cdm1 <- glm(sat ~ TreatGroup, data = d.noext, family = binomial(link = "logit")); summary(h1cdm1)
h1cdm2 <- glm(sat ~ TreatGroup + BG2_Age + BG2_Sex_1 + Dem3_Edu_lev + Dem3_Occpn_if_desk + BG2_City + Knowl9_IMAGINE + act_qrt + AddOut8_KnowCh + 
                PreT4_Context_PoliKnowle, data = d.noext, family = binomial(link = "logit")); summary(h1cdm2)
stargazer(h1cdm1, h1cdm2, single.row = T, digits = 2, star.cutoffs = c(.05, .01, .001,.0001))

########################################
########################################
##Table A9: Evaluations of INGO and SOE Effectiveness
########################################
########################################

#Run regressions reported in Table A9
reg_pool_3 <- lm(PosTre6_Brand_Eff_EuroNGO ~ brand_pooled + ext_pooled + brand_pooled*ext_pooled,data=d)

reg_pool_3c <- lm(PosTre6_Brand_Eff_EuroNGO ~ brand_pooled + ext_pooled + brand_pooled*ext_pooled + PreT4_Eff_INGO +
                    BG2_Age + BG2_Sex_1 + Dem3_Edu_if_col + Dem3_Occpn_if_desk + BG2_City + Knowl9_IMAGINE + act_qrt + 
                    AddOut8_KnowCh + PreT4_Context_PoliKnowle,data=d)

reg_pool_4 <- lm(PosTre6_Brand_Eff_ChinaSOE ~ brand_pooled + ext_pooled + brand_pooled*ext_pooled,data=d)

reg_pool_4c <- lm(PosTre6_Brand_Eff_ChinaSOE ~ brand_pooled + ext_pooled + brand_pooled*ext_pooled + PreT4_Eff_ChinaSOE +
                    BG2_Age + BG2_Sex_1 + Dem3_Edu_if_col + Dem3_Occpn_if_desk + BG2_City + Knowl9_IMAGINE + act_qrt + 
                    AddOut8_KnowCh + PreT4_Context_PoliKnowle,data=d)

stargazer(reg_pool_3, reg_pool_3c, reg_pool_4, reg_pool_4c, single.row = T, digits = 2, star.cutoffs = c(.05, .01, .001,.0001))

########################################
########################################
##Table A10: Attribution, Non-Pooled Sample
########################################
########################################

#Run regressions reported in Table A10
h5m1a <- glm(ingo_eff ~ TreatGroup, data = dext, family = binomial(link = "logit"))

h5m2a <- glm(ingo_eff ~ TreatGroup + BG2_Age + BG2_Sex_1 + Dem3_Edu_lev + Dem3_Occpn_if_desk + BG2_City + Knowl9_IMAGINE + act_qrt + AddOut8_KnowCh + 
               PreT4_Context_PoliKnowle, data = dext, family = binomial(link = "logit"))

h5m1b <- glm(soe_eff ~ TreatGroup, data = dext, family = binomial(link = "logit"))

h5m2b <- glm(soe_eff ~ TreatGroup + BG2_Age + BG2_Sex_1 + Dem3_Edu_lev + Dem3_Occpn_if_desk + BG2_City + Knowl9_IMAGINE + act_qrt + AddOut8_KnowCh + 
               PreT4_Context_PoliKnowle, data = dext, family = binomial(link = "logit"))

stargazer(h5m1a, h5m2a, h5m1b, h5m2b, single.row = T, digits = 2, star.cutoffs = c(.05, .01, .001,.0001))

########################################
########################################
##Figure A5: Effect of Moving from Single-Actor to Multi-Actor on Project Satisfaction (Figure A8 in conditionally accepted version of manuscript)
########################################
########################################

#Create effects plot for Figure A5
Externality <- c("Positive","Positive","None","None","Negative","Negative")
Branding <- c("INGO","SOE","INGO","SOE","INGO","SOE")
ATE <- c(.1145,.3049,.3476,.3,.0735,.059)
LCI <- c(-.1781,0.0061,0.06,.01486,-0.2168,-0.248)
HCI <- c(.4072,.6038,.6357,.5851,.3638,.366)
d.ate.plot <- as.data.frame(cbind(Externality,Branding,ATE,LCI,HCI))
d.ate.plot$ATE <- as.numeric(as.character(d.ate.plot$ATE))
d.ate.plot$LCI <- as.numeric(as.character(d.ate.plot$LCI))
d.ate.plot$HCI <- as.numeric(as.character(d.ate.plot$HCI))

f.a5 <- ggplot(d.ate.plot, aes(x=Externality,y=ATE, color=Branding)) +
  geom_point(position=position_dodge(width = .5),size=5,aes(shape=Branding)) +
  geom_errorbar(aes(ymin=LCI, ymax=HCI,color=Branding),position=position_dodge(width = .5),width=.2) +
  ylim(-.3,.85) +
  coord_flip() +
  geom_hline(yintercept=0, linetype="dashed", color = "black") +
  scale_color_manual(values=c("black", "black")) +
  theme(axis.title.x = element_text(size=16),
        axis.text.x  = element_text(size=14,color="black"),
        axis.text.y  = element_text(size=14,color="black"),
        axis.title.y = element_text(size=16),
        legend.title=element_text(size=16),
        legend.text=element_text(size=14),
        legend.key.size =unit(2,'cm')) +
  ylab("Average Treatment Effect Estimate (Raw 1-7 scale)")

##Close .log file
closeAllConnections()
